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1.  Introduction 

In  recent  years,  widespread  attention  has  been  given  to  the  computation  of 
stationary  distributions  of  Markov  chains.  A  variety  of  methods  have  been 
suggested  and  implemented.  Before  considering  an  alternative  way  for  finding 
such  distributions,  it  is  of  interest  to  give  a  brief  survey  of  the  techniques 
that  have  been  employed. 

Paige,  Styan  and  Wachter  (1975)  presented  a  comprehensive  survey  of  eight 
different  algorithms  involving  a  variety  of  procedures  including  the  use  of 
generalized  inverses,  rank  reduction,  least  squares  and  power  methods.  Their 
recommendation  was  a  direct  method  that  involved  transforming  the  singular  set 
of  stationary  equations  into  a  non-singular  system  using  a  rank  one  modification 
followed  by  Gaussian  elimination  with  row  pivoting.  A  further  study  by  Harrod 
and  Plemmons  (1984)  provided  another  direct  approach  based  upon  the  LU 
factorization  using  Gaussian  elimination  without  pivoting. 

Iterative  techniques  and  approximation  methods  have  been  surveyed  by  Koury, 
McAllister  and  Stewart  (1984).  When  the  transition  matrix  is  large  and  exhibits 
a  nearly  completely  decomposable  structure  it  is  shown  that  a  method  of 
"aggregation"  can  be  combined  with  point  and  block  iterative  techniques  to 
produce  methods  which  converge  rapidly  to  the  stationary  probability  vector. 

Sheskin  (1985)  presented  a  partitioning  algorithm  that  used  a  matrix 
reduction  routine  that  partitions  the  transition  matrix  to  create  a  sequence  of 
smaller  order  transition  matrices  followed  by  a  vector  enlargement  routine  that 
enables  the  components  of  the  steady  state  vector  to  be  determined  sequentially. 
A  related  procedure  was  developed  by  Grassmann,  Taksar  and  Heyman  (1985)  using 
the  theory  of  regenerative  processes.  They  derived  relationships  between  the 
steady  state  probabilities  which  are  then  used  to  develop  a  numerical  algorithm 
to  find  these  probabilities.  Both  of  these  latter  two  techniques  appear  to  be, 
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in  effect,  modifications  of  Gaussian  elimination. 


More  recently,  Meyer  (1987),  has  utilised  the  concept  of  "stochastic 
complementation"  whereby  an  irreducible  chain  is  uncoupled  into  smaller 
irreducible  chains  whose  stationary  distributions  can  be  coupled  back  together 
to  produce  the  stationary  distribution  of  the  original  chain. 

In  this  paper  an  entirely  new  approach  involving  the  analysis  of  perturbed 
Markov  chains  is  considered.  In  Hunter  (1986)  techniques  for  updating  the 
stationary  distribution  of  a  finite  irreducible  Markov  chain,  following  a  rank 
one  perturbation  of  its  transition  matrix,  were  presented.  In  this  current 
paper,  these  techniques  are  utilised,  to  construct  a  general  procedure  for 
determining  the  stationary  distribution  of  any  finite  irreducible  Markov  chain. 

A  significant  feature  of  the  proposed  algorithm  is  that  at  no  stage  does  a 
system  of  linear  equations  have  to  be  solved  and  consequently  there  is  no 
reliance  upon  computer  subroutines  for  matrix  inversion  or,  the  more  generally 
accepted  method  of  solution,  Gaussian  elimination  with  or  without  pivoting. 

The  basic  idea  is  very  simple.  Suppose  the  steady  state  probability  vector 
it'  of  an  m-state  irreducible  Markov  chain  with  given  transition  matrix  P  is 

required.  Let  P  be  the  transition  matrix  of  smother  irreducible,  m-state, 
o 

Markov  chain  with  known  stationary  probability  vector  ir'.  By  replacing, 

~v) 

successively,  the  elements  of  each  row  of  PQ  with  the  corresponding  row  elements 

as  specified  by  P  and  recomputing  the  stationary  probability  vector  of  the 

resultant  perturbed  transition  matrix,  the  vector  tt’  can  be  transformed,  in  m 

~o 

stages,  to  tt' ,  by  a  series  of  m  updates. 

As  the  irreducibi 1 i ty  of  a  Markov  chain  is  governed  by  the  location  of  the 

positive  entries  in  its  transition  matrix,  to  ensure  the  i rreducibi 1 i ty  of  each 

perturbed  Markov  chain  it  is  sufficient  to  commence  with  P  containing  positive 

o 

elements  placed  at  least  in  the  same  position  as  those  in  P. 


Consider  starting  with  the  trivial  doubly  stochastic  matrix  Pq  with  each 
element  having  the  value  1/m,  so  that  P  =  ee'/m,  where  e'  =  (1,1 . 1)  is  a 

O  »v»v  'V 

vector  of  ones.  As  can  be  easily  shown,  v'  =  e'/m. 

'vO  <v 

For  i=l,2,...,m,  let  e^  be  the  i ^elementary  (column)  vector  with  a  one  in 
the  i t^1  position  and  zeroes  elsewhere.  Let  =  e^P  be  the  i t^1  row  of  P  and  let 


P.  =  P.  ,  +  e.b;, 
1  l-l  ^1^,1 


where 


b ;  =  p ;  -  e'/m. 

/■v,!  Zi  ^ 


Let  7rj  be  the  stationary  probability  vector 


chain  with  transition  matrix  P  ,  and, 
required  vector  ir' . 


since  P  = 
m 


(1.1) 


associated  with  the  Markov 
m 

2  e.p!  =  P,  ir'  is  in  fact  the 


2.  General  Theory 

The  construction  of  the  algorithm  is  based  upon  the  following  results. 
Theorem  2.1=  Let  P^  be  the  transition  matrix  of  a  finite  irreducible  Markov 
chain  with  stationary  probability  vector  ir!  . 

>vl 

(a)  I  -  P.  +  t.ul  is  non-singular  if  and  only  if  u|e  ?  0  and  ir! t.  ^  0. 

1  1  /V/  ^  1  /V  /vIa/I 

(b)  Under  the  conditions  of  (a), 

jr'  =  a'/a'e,  (2.  1) 

r\j  ^  1  (v  *  <v 

where  a'=u'[I-P+t.u!]*.  (2.2) 

Proof :  For  (a)  see  Theorem  3.3  in  Hunter  (1982)  and  for  (b)  see  Corollary  4.1.2 
in  Hunter  (1982).  D 


>  .  • . 


Theorem  2.2  If  X  is  non-singular  and  b‘X  a  ^  -1,  then 


X  *ab  X  1 


(X  +  ab' )  1  =  X  1 


1  +  b‘X  *a 


(2.3) 


Proof .  This  is  the  Sherman-Morr ison  formula.  See  Golub  and  Vein  Loan  (1983)  p3. 


Suppose  that,  following  the  i  perturbation,  the  stationary  probability 


vector  ir!  has  been  found  for  the  Markov  chain  with  transition  matrix  P. ,  as 

l 


given  by  (1.1),  by  using  the  procedure  described  by  Theorem  2.1(b)  for  suitable 


choices  of  t.  and  u.. 

^.i  ^,1 


In  Hunter  (1986)  it  was  shown  that  it  is  possible  to  find  an  expression  for 


it'  . ,  associated  with  P ...  using  the  same  procedure  outlined  in  Theorem  2.1(b), 
1*1  l  *  X 


by  choosing  the  t^+^  and  u^+^  in  such  a  way  that  [I  -  P.+1  +  ti+iui+i^  1  can 


determined  from  the  earlier  deduced  expression  [I  -  P  +  t  u']  ,  without 

1  (V  •  X 


performing  an  additional  matrix  inversion. 


For  the  particular  situation  under  consideration,  for  i=0,l . m-1,  if 


t, , .  =  e.,  and  u.,,  =  u.  +  b.  ,,  where  b.  is  given  by  (1.2),  then,  from  (1.1), 

^.i+1  ^i+1  ^i+l  ^i  ^.i+1  ^i  &  v  ’  v  ' 


I  -  P...  +  t .  u'  .  =  I  -  P.  +  e.^.u!, 
i+1  ^i+l^i+1  i  ~i+l^i 


-  1  -  pi  * 


(2.4) 


if  [I  -  Pj  +  t^ul]  exists,  from  the  proof  of  Theorem  3.3  in  Hunter  (1982), 


ui[I  ‘  pi  +  Vi^  1 


(2.5) 


BS 


J*  ^  "V" 


’a  ^  W~*T  W-^K  —  ■  r\  *.  - 


Thus,  using  (2.3).  (2.4)  and  (2.5),  if  A.  =  [I  -  P.  +  t.uj]  exists. 


A.  ,  =  A.  [I  +  ( t .  -  t .  .)  - ; - ] 

1  +  1  1L  v^,l  ^l+ly  w.  t.  ,  J 

~1^.1+1 


(2.6) 


Equation  (2.6)  is  ideally  suited  for  recursive  operations  once  an  initial 

inverse  A  =  [I  -  P  -  t  u’]  *  has  been  determined.  However,  because  of  the 
o  o  ^,0^,0 

form  of  P  that  has  been  selected,  if  t  =  e  and  u'  =  e'/m,  no  matrix  inverse 

o  ^o  ^  ~,o  ^ 

has  to  be  computed  since,  in  this  instance. 


I-P  +  t  u'  =  I  -  ee'/m  +  ee'/m  =  I. 
o  ~o^o  —  _ 


Furthermore,  using  (2.2), 


and,  f rom  (2.1), 


a '  =  u '  =  e ‘ /m , 

/N/O  ^vO  >v 


7 r '  =  e ' /m . 

~o  ^ 


The  basic  algorithmic  procedure  now  follows-' 

Let  t  =e,  u  ’  =  e'/m,  A  =  I,  ir"  =  e'/m. 

A/O  "V  /vO  -V  O  /V. 

For  i=l,2 . m,  let  t.  =  e.  ,  and  u’.  =  u!  +  p!  -  e\ 

'V.  ^  <vl  <vl  1  /%* 


Compute 


Compute 


Compute 


i.  =  A  ,[I  +  (t.  ,  -  t.)ir;  ,1.1. 

i  i-lL  \,i-l  ^l'^i-l  ^i-l^i J 


a ;  =  u: a. . 

^i  i 


ir  :  =  a !  /a !  e . 


(2.7) 


(2.8) 


(2.9) 


Then  ir'  =  ir'  is  the  stationary  probability  vector  of  the  Markov  chain  with 
transition  matrix  P  =  [p  ]. 

Since  the  elements  of  any  stationary  probability  vector  are  always 
positive,  Trlt  and  ir'  t.  are  both  positive.  Further,  by  induction,  for 

<vl  —  I«vi 

i=0 ,  1 . m , 


u !  e  =  1 , 


(2.10) 


so  that  the  conditions  of  Theorem  2.1  and  2.2  are  satisfied. 
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3.  Refinements  to  the  algorithm 

Although  the  procedure  suggested  in  Section  2  will  lead  to  the  required 
stationary  probability  vector  there  are  some  modifications  that,  if  employed, 
will  lead  to  a  more  efficient  procedure. 


3.1  Modification  to  the  ir!  computation 
-  - - 

The  ultimate  aim  of  the  algorithm  is  to  determine  ir*  =  ir' .  Unless  the 

~  ~m 

stationary  distributions  of  the  intermediate  perturbed  Markov  chains  are 
required,  some  simplification  can  be  effected  by  observing  that  (2.7)  requires 
ir'  through  its  scaled  version  ir'  /ir! t .  .  The  scaling  suggested  by  (2.9)  is 

not  required  until  the  final  step  when  i=m. 

Thus,  for  i=0, 1 , . . . ,m-l ,  let 


u!  =  ir ! /ir !  t .  , ,  . 
^.i  ~i~i+l 


(3.1) 


Then,  o'  =  e'  and  for  i=l,2 . m-1, 

A/O  (V 

u!  sulA./ulA.e.,,  . 

~i  i  i^i+l 

If  ir!  is  required  then  it  can  be  recovered  simply.  For  i=0, 1 . m-1 

A/  1 


(3.2) 


ir;  =  v'./v'.e 


(3-3) 


At  the  final  step,  compute  ir^  using  (2.8)  and  (2.9). 


3.2  Modification  to  the  A^  computation 

With  the  notation  introduced  in  Section  3.1,  it  is  easily  seen  that  the 
early  terms  in  the  (A^)  sequence  are  given,  after  simplification,  as 


A  =  I. 
o 

A.  =  A  +  (e  -  e  )u' 
1  o  ~l~o 


A2  "  A1  +  (®1  “  ^2^1  ' 


A3  "  A2  +  ^  U13^®1  +  P13!!2  ’ 


(3.4) 


(3-5) 


(3.6) 


(3.7) 


A4  ”  A3  +  u14  ^  u13^u24^1 

+  ^14  '  U13U24}!2  +  U24^3  *  !4]S 

where  v.  .  =  u!e.,  so  that  v .  .  .  =  1  for  i=l,2 . m-1. 

ij  ~i~J  1,1+1 

The  above  results  provide  motivation  for  the  following  theorem. 


Theorem  3 . 1 •  For  n=0,l . m-1 


where 


A  ,  =  A  +  B  , 
n+1  n  n 

B  =  A  ( e  -  e  , ) u '  =  b  i>' 
n  nvn  ^n+l  Toi  ~n  ~n 


with  e  =  e,  so  that  b  =  e  -  e  and  for  n=l,2 . m-1, 

^o  ^  ^.o  ~1 

b  =  b.  e,  +  ...  +  b  e  -  e  ,  . 
liul  mun  ~n+l 


(3.8) 

(3.9) 

(3.10) 


Proof :  The  theorem  is  obviously  true  by  inspection,  from  (3.4)  to  (3.7),  for 

n=0, 1.2,3.  Assume  that  (3.9)  and  (3.10)  hold  for  n=0,l,...,k  so  that 


Hence 


\+!  -  \  +  \  =  Vi  +  Vi +  V 

k 

=  ...  =  I  +  2  B  . 

~  n 
n=0 


^+1  '  "V+i^+l  5t+2}S+l 


(3.11) 


implying  that 


-  (I  +  Bq  +  2  Bn)(ek+,  ^+2)^+1 

n=l 


k  n 


5k+i  =  ci  *  (« -  J,bmns  '  '  W 

n=l  m=l 

Since  !'(!k+i  "  !k+2}  =  0i 


k  n 


Vl  =  !k+l  "  !k+2  +  2  (  2  bmn!m  “  !n+l){Vk+l  '  Vk+2 

n= 1  m= 1 
k  k 

=  ^k+1  ~  tk+2  +  2 /  2  bmn^Un.k+l  “  un,k+2^fm 

m= 1  n=m 


)• 


k+1 
2  \ 
m=2 


2  ^m-l.k+1  Vl.k+2^  ' 


showing  that  is  a  linear  combination  of  . ^k+2  wit^  ^k+2  ^iaving 

coefficient  -1.  Thus  (3.10)  is  true  for  n=k+l  and  the  theorem  follows  by 
induction. 


Observe  that 


b.  ,  =  2  b  ,  ,  e  -  e,  _ , 

Tk+1  ,  m,k+l^m  ^k+2 

m=  1 


where 


and  for  m=2 . k, 


b  1 , k+ 1  2  bl.n^n.k+l  L’n,k+2^  1 

n=l 

bk+l.k+l  =  Uk. k+2  '  (since  ^k.k+r1^ 


a  ,  . ,  =  2  b  (d  .  ,  -  v  .  , ) 
m,k+l  .  m,nv  n,k+l  n,k+l' 

n=m-l 


(3. 13 
(3.  h; 

(3.15)  c 


Note  that 


B  =  b  v' 
n  „ji  ^n 


In  [u  . . v 

L  nl  n2 


,  u  ]  , 

nmJ 


’ln^nl 

blnUn2 

b  v 

In  nm 

*2nUnl 

b2nUn2  ' ' 

b„  u 

2n  nm 

nnUnl 

bnnUn2 

b  u 
nn  nm 

_Unl 

0 

_Un2  ' ’ 
0 

-u 

nm 

0 

0 

0 

0 

(3- 10) 


Thus,  in  the  matrix  B  ,  all  the  entries  in  rows  numbered  n+2 . m  are  zero. 

n 

Obviously,  this  has  considerable  significance  in  the  calculation  of  the  matrices 
A.  (i=2,...,m)  as  required  in  the  algorithm. 


wmmmmmmmmmmm 


The  updating  process,  given  by  (2.7),  can  be  replaced  for  i=2 . m,  by 


A .  =  A .  ,  +  B .  ,  , 
l  i-l  l-l 


(3. 17) 


where 


B.  .  =A.  ,(e.  ,  -  e .  )i> !  , 
i-l  i-l^i-l  ^i'^i-I 


(3. IS) 


is  such  that  only  the  first  i  rows  of  B  ^  are  computed  with  the  remaining 


entries  set  equal  to  zero. 


Furthermore,  from  (3.16),  some  of  the  rows  of  B^  have  a  special  form  and  do 


not  require  computation.  In  particular  the  (n+l)th  row  is  simply  -v '  while. 


from  (3.14),  the  n  row  is  v  ,  ,  times  v‘ . 

n- 1 , n+ 1  ^n 


Note  also  that  the  (n+l)st  column  of  B  is  b  since  v  ,=] 

n  ^n  n ,  n+ 1 


There  are  also  some  other  checks  that  can  be  applied. 


Theorem  3.2:  For  i=l,2 . m, 


A .  e  .  =  e , 
i~i  ~ 


(3.19) 


1  'A 

—  e  A.  =  e, 

rn  ^  i  r\. 


(3.20) 


e'B.  =  O' . 
^  i 


(3.21) 


Proof :  Since  A^  =  [I  -  +  t^uj]  .  equation  (3.17)  of  Hunter  (1982)  implies 


A .  t .  =  — ; =  e  , 

i~i  u.e  ^ 

~i~ 


yielding  (3.19)  with  t.  =  e^. 


Equation  (3.20)  is  obviously  true  when  i=l  since 


ZK\  =  !'[I  +  (!  ‘  !i}!’] 


=  e ‘  +  (m-1 )e ’  , 


=  m  e  . 


Thus,  by  induction,  if  (3.20)  is  true  for  i=l,2 . k,  from  (2.7)  and 


(3.1). 


V.W  V.'/.VV.V.V  ' 


-V  *  -  »  m*V  m  ft  - 


•/--.v.v.v.v. 


m  Z  K+l  m  ~  TcL  .Tk+l'3c- 


“  +  <3t  "  Sc+i>ik]  • 


=  e  '  +  ( e  1  e,  -  e  '  e,  ,  )  v,‘  =  e ' 
~  v _ k _ k+lyJ( 


Thus  (3.20)  and  hence  also  (3.21)  follow. 


A  consequence  of  Theorem  3.2  is  that  the  i  column  of  consists  solely 

of  unit  elements  while  the  sum  of  the  elements  of  each  column  of  B.  is  zero. 

i 

3.3  Modification  to  the  a.  computation 

Although  the  {A^}  (i=0,l . m)  sequence  plays  an  integral  role  in  the 

procedure,  the  matrices  A^  are  required  only  to  obtain  the  sequence  of  vectors 
aj  =  u^A^  and  hence  the  vectors  ul.  Thus  is  is  worth  examining  whether  it  is 

possible  to  dispense  with  explicit  calculation  of  the  A^  by  deriving  the  (a^) 


(i=l,2 . m)  sequence  recursively. 


Theorem  3.3:  For  i=0, 1 ,2, . . . ,m-l , 


a'  =  i>!  -  e '  +P’,iA.l1 
^l+l  ^i  ^  Ti+1  l+l 


Proof :  First  observe  that,  from  (2.8), 


;i  -  “iAi  -  EiAr 


so  that  (3.22)  holds  for  i=0  since  uq  =  e ' . 

In  general,  for  i  =  1 , 2 . m-1.  from  (2.8), 

;Ui  -  “h iAn-r 


“  [uj  -  (e'/m)  +  p;tl]Aitl 

“iVi  =  “TAi  *  Ai(si  -  ti+i^d- 
=  ;i  *  Zi<Zt  ~  ti+i’h' 


=  a',  +  (“n  -  a,  ,,,)«:  . 


where  a.  .  =  a'.e..  But,  from  (3.19)  and  (2.10), 
ij 

a.  .  =  u ! A . e .  =  u ! e  =  1  , 
n  ^.i  i^i  ^,1^. 

and,  since  from  (3.2)  and  (2.8),  uj  =  aYa.  ^  +  (3-24)  becomes 


(3.25) 


uiA.  ,  =  i>; . 
~i  i+i  ^i 


(3.26) 


Equation  (3.22)  now  follows  from  (3.23)  upon  substitution  of  (3.26)  and 


(3.20). 


Theorem  3.3  shows  that  in  updating  from  a',  to  a'  the  term  p!,.A.  must 

*  °  ^,1  ^i+l  Ti+1  i  +  l 

be  computed.  The  calculations  of  a(  ,a' . a!  require,  successively,  p.' . p! 

and  for  a!  .  this  is  the  first  time  p!  ,,  the  (i+l)th  row  of  P  is  involved. 

~i+l  Ti+1 

Although,  pj+^A^+^  can  be  expressed  in  terms  of  A^ ,  as  can  be  seen  from  the 

next  theorem,  very  little  advantage  is  gained  since  such  terms  are  required  for 
each  i=0,  1 . m-1 . 


Theorem  3.4:  For  i=0, 1 . m-1 


p!  ,A.  ,  =  p'  A.  +  i>!  -  (p !  ,  . A .e .  .  ,  )i> ! 

£i+l  i+l  iCi+1  i  ~i  v£i  +  l  i^i+lO 


(3.27) 


Proof :  For  i=0. 


p1Al  =  P^1  +  (e  "  ei)e'l  • 

=  Pi'  +  (Pie)e  "  (Piei  )e 

/V*  /v^rv  >V  <V  1  »V  A  ^ 


and  the  result  follows,  since  p.'e  =  1.  v '  =  e'  and  A  =  I. 

^o  ^  o 

In  general,  for  i=l,2 . m-1,  from  (2.7)  and  (3.1), 

Ei+iAi*i  =  h+i[Ai  +  V't  ■ 

Equation  (3.27)  follows  since,  from  (3.19), 

p'A.e.  =  p!  •  =  1 . 
r.i+1  i~i  ~i+i~ 


As  a  consequence  of  Theorem  3.3,  it  is  suggested  that  (2.8)  in  the 
algorithm  be  replaced  by  (3.22). 

3.4  Modification  to  the  w'  =  ir‘  computation 


At  the  final  step  of  the  algorithm  Am  can  be  computed  and  consequently 

derived  as  a'/a’e  where  a’  =  u‘A  .  However,  A  need  not  be  explicitly 
~jn  ^rru.  ~m  „jn  m  m 

determined  since,  from  Theorem  3.3  and  3.4, 


a'  =  v'  .  -  e '  +  p '  A  , 
^jn  „jn-l  ji  m 


sre  p '  A  =  p  '  A  ,+u’ 

2m  m  2m  m-1  ^jn-1 


(p ' A  ,  e  )v  '  ,  . 
v2m  m-l^jn'-jn-l 


4.  Recommended  procedure 

As  a  consequence  of  the  refinements  discussed  in  Section  3,  it  is  suggested 
that  the  algorithm  be  constructed  as  follows: 


1. 

Let 

Aj  =  I  + 

2. 

Let 

;i  =  hA 

3. 

For 

i=1.2. . . 

(a)  v !  =  a'./a’.e .  .  , 

<b>  Bi  -  -  It.dii  • 

(c)  A1+1  =  a.  ♦  B,  . 

(d)  ei.i  =  ;i  - +  eLiaui 


4 .  Let  u '  ,  =  a ’  /a '  , e 


5.  Leta'=2i;'-e'+p,A  ,  -  (p'A  ,e  )u’  . 

^m-1  ^  2m  m-1  2m  m-l^m^^m-1 


6.  Let7r'  =  ir '  =  a ' /a ' e  . 

A/  aJH  aJIXv 


vY#  V  V-V'/v.  /V.  •.*  ,  /a  mi'.  \itohr*  >>  hNAftaN 


The  order  of  the  number  of  arithmetic  operations  (multiplication  and 
division)  required  to  determine  ir'  can  be  estimated  as  follows.  The  computation 


of  the  B.  and  the  P-A.  have  a  dominant  effect  on  the  number  of  operations 

required.  Since  A.(e.  -  e.  ,)  is  effectively  the  difference  of  two  columns  of 

A.,  only  mi  operations  are  required  to  determine  B^ ,  taking  into  consideration 
that  Bj  has  only  (i+1)  non-zero  rows,  and,  as  a  consequence  of  (3.21),  that  the 
elements  of  one  row  can  be  found  from  the  other  rows  using  the  fact  that  each 
column  sums  to  zero.  On  the  other  hand,  for  a  general  transition  matrix,  p^A. 

2 

will  require  m  operations,  although  this  can  be  reduced  to  m(m-l)  since  the  i th 
element  of  this  row  vector,  p'.A.e.  =  p.'e  =  1,  (by  using  (3.19)).  Since  the 

other  calculations  required  are  relatively  insignificant  in  comparison,  the 
total  number  of  operations  is  of  the  order  of 

m- 1  m  ~  ^ 

I  mi  +  2  m(m-l)  =  3m  (m-l)/2.  i.e.  of  order  3m  /2. 
i=l  i=l 

To  solve  for  the  stationary  distribution  directly  using  Gaussian 

3 

elimination  requires  of  the  order  of  4m  /3,  while  to  solve  directly  using  a 

3 

matrix  inversion  routine  requires  of  the  order  of  2m  operations,  (see  Isaacson 
and  Keller  (1966)). 

The  procedure  is,  in  effect,  finding  the  stationary  distribution  of  m 
different  irreducible  Markov  chains  and  consequently  the  routine  that  has  been 
developed  offers  much  more  information  than  other  techniques  currently 
avai lable. 

Although  it  has  been  suggested  that  the  algorithm  proceed  row  by  row,  there 
is  no  necessity  to  adhere  to  a  strict  sequential  ordering  of  the  rows.  The 
procedure  as  outlined  by  (2.7),  (2.8)  and  (2.9)  can  easily  be  adapted  to  such 
changes  by  altering  the  t^  and  uj.  A  consequence  of  this  is  that  the  effect  of 


4 


changing  selected  transition  probabilities  upon  the  stationary  distribution  can 
easily  be  determined.  (See  also  Hunter  (19S6)). 

The  procedure  also  offers  the  opportunity  to  utilise  the  structure  of 
special  transition  matrices.  For  example,  if  the  transition  matrix  of  the  chain 


is  banded  with  p . ^  =  0  for  j<i-g  and  j>i+h.  which  occurs  in  some  queueing 


models,  the  calculation  of  will  require  at  most  (g+h)m  operations  and  the 


3  2 

algorithm  will  require  on  the  order  of  only  m  /2  +  (g+h)m  operations. 


5.  Structural  results 

In  Section  3.2  expressions  for  the  first  few  terms  of  the  (A^  sequence 
were  derived.  By  using  those  terms  and  working  through  the  first  few  steps  of 
the  algorithm  it  can  be  shown,  that,  following  simplification,  for  i=l,2,3, 


a ;  = 

^1 

(p.  e' 
v  io^ 

+  PilPi 

+  . 

(5.1) 

"i  = 

(P-  e ' 

1  io^ 

+  PilPi 

+  . 

*  “itEp^lo  *  uilpl 

,  i+1  + 

•  +  pii> 

Pi.  i+P 

.  (5.2) 

II 

-  *r-4 

►=  ? 

+  pnpi 

+  . 

*  uiiEi)x(mMio  *  Mii 

+  .  .  .  + 

Pit)  1 

(5.3) 

where 

p10  = 

1  •  pir 

pn  = 

l. 

p20  = 

(1  '  pn 

)(1 

- 

p22*  "  p12p2r 

P21  = 

1  +  P21 

-  p 

22  ‘ 

U22  = 

1  "  P11 

+  p 

12 

- 

p30  = 

(1  -  pn 

)(1 

- 

P22^ 1  p33^  P12P23P31 

P13P21P32 

- 

p13^ 1  ’ 

P22 

)p31  ^  P11^P23P32  + 

P 1 2P2 1 

( 1  ~  p33^ 

• 

P31  = 

P21^  - 

P33 

+ 

p32^  +  ^ 1  p22^1 

P33  +  P 

31 ^  +  P23 

(p31  * 

P32^ 1 

P32  = 

P3 1 ^ P 1 2 

-  P 

13^ 

'  +  p32(1  “  P11  +  P13> 

+  (1  - 

P11  + 

P12^  ’ 

p33  = 

0  -  pn 

)(1 

- 

p22  +  p23^  +  P12^P23 

-  p21) 

+  P13{1  " 

p22  + 

p21}- 
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The  general  structure  exhibited  by  (5.1),  and  hence  also  by  (5.2)  and 

(5.3),  holds  for  all  i=1.2 . m.  [A  proof  by  induction  shows  that  if  (5.1)  and 

(5.2)  hold  for  i=l . n  then,  since 

p'A  ,  =  p '  [I  +  B  +...+B], 

Cn+1  n+1  Cn+1 L  o  nJ 

=p'  , [I  +  b  u ‘  +  ...  +  b  u  '  ]  , 
iji+1 L  ^o^o 

=  p’  ,  +  ( p ‘  ,b  )u'  +  ...  +  (p'  .b  )u'  , 

Cn+1  vCn+l~o'^o  vOn+l~nOn 

using  (3.22)  with  i=n,  a'  is  a  linear  combination  of  v ' . u ' ,  p'  ,  i.e.  of 

v  ’  ^n+1  Cn+1 

e',p’,...,p'  Furthermore,  the  coefficient  of  p'  ,  is  unity  whereby 

~  Cl  Cn+1  Cn+1 

establishing  the  general  structure  of  oh^.] 

Note  also  that,  from  (5.1)  and  (3.25),  ahe^  =  1,  for  i=l . m  and  thus 


Pio  +  PilPli  +  ••  +  PiiPii  =  Pii  • 


Further,  for  i=l,2  it  can  be  shown,  by  direct  verif ication,  that 


Pio  +  PilPl,i+l  +  +  PiiPi,i+l  =  Pi+l,i+l 


(5.4) 


which  implies  that 


Uii  =  ”i*i  =  Pii/Pi+I,i+1 


a.  ...  -a !e,  ,  =  u.  ,  .  ./p. . 
i,i  +  l  ^i^i+1  l  +  l,  l+l  li 


(5.5) 


(5.6) 


(5.7) 


results  that  it  has  not  been  possible  to  establish  in  general. 

Let  (I-P)j  be  the  leading  i th  order  principal  submatrix  of  I-P  formed  by 


deleting  all  but  the  first  i  rows  and  columns  then,  for  i=l,2,3. 


pio  =  det(I-P)1 


(5.8) 


« 


For  the  special  case  when  m=3,  with  the  notation  used  earlier  in  this 


section,  it  can  be  verified  that 


a)  =  (4. , .  p.n,  p._)/p. .  . 
^,1  1  1 1  i2  i3'  11 


"i  =  (pir  4i2-  4i3^/Mi+i,i+i 


!i  =  {4il'  Mi2’  pi3)/(pil  +  pi2  *  pi3>' 


(i=l ,2,3) , 

(i=1.2), 

(1=1. 2. 3). 


(5.9) 

(5.10) 

(5.11) 


where 


412  ”  1  P11  +  P12  "  P22‘ 


p13  =  1  -  pll  +  p13 ' 


p23  “  p33  "  ^30 ' 


Observe  that  tt'  it'  and  ir'  give,  respectively,  the  stationary  probability 

A/  A 

vectors  of  the  Markov  chains  whose  transition  matrices  are 


P11  p12  P13 


P21  p22  p23 

P31  p32  p33 


In  examining  (5.11),  with  i=3,  it  can  be  shown  that  4.  .  =  3D.,  (j=l,2,3). 

3  J  J 

where  is  the  determinant  formed  by  striking  out  the  jth  row  and  j th  column  of 
I-P.  This  leads  to  an  expression  for  the  stationary  probability  vector  of  a 
general  irreducible,  three  state,  Markov  chain  as 


pll 

p12 

P13 

’  P11 

P12 

p13 

1/3 

1/3 

1/3 

P21 

P22 

P23 

1/3 

1/3 

1/3 

1/3 

1/3 

1/3 

W  =  (D  D  D  )/  2  D  . 

1  Z  j=l  J 


(5.12) 


The  natural  extension  of  (5.12)  for  a  general  finite  irreducible  Markov 
chain  is  also  true,  such  a  result  being  attributed  to  Mihoc  by  Frechet  (1950) 
and  rediscovered  by  Singer  (1964). 
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Although  the  full  details  of  a  proof  of  the  generalization  of  (5.12)  using 

the  techniques  of  this  paper  have  not  been  worked  out,  it  is  conjectured  that 

for  an  m-state  chain  p  .  =  mD . .  (a  result  that  holds  for  m=2,3),  so  that  the 

J 

procedures  proposed  in  this  paper  appear  to  lead  to  an  effective  algorithmic 
construction  of  Mihoc's  technique. 

6,  Final  comments 

The  initial  choice  of  P  as  ee'/m  ensures  that  it  is  possible  to  start  with 

o  —  r 

an  irreducible  Markov  chain  whose  stationary  distribution  is  easily  found 
without  having  to  compute  a  matrix  inverse  or  to  solve  a  general  set  of  linear 
equations.  The  fact  that  every  element  of  Pq  is  specified  leads  to  a  sequence 
of  matrices  that  are  "dense".  Is  it  possible  to  start  with  a 

different  Markov  chain,  say  one  that  is  relatively  sparse,  whose  stationary 
distribution  is  well  known  and  such  that,  for  the  early  recursions,  the 
equivalent  sequence  A^.Ag....  retains  such  a  sparsity  property? 

The  periodic  Markov  chain  with  entries  p(?\=l,  (i  =  l,2 . m-1).  and  p^°^=l 

is  a  potential  candidate  for  Pq,  whose  stationary  probability  vector  is  also 
tt'  =  e'/m.  Even  if  t  and  u  can  be  specified  so  that  A  =  f I  —  P  +  t  u'l  1 

^  O  O  ^.0^,0 

has  a  simple  structure  much  care  would  be  required  in  carrying  out  any 

sequential  row  modification  with  this  P  .  For  example  if,  for  the  specified  P 

transition  matrix,  p^  =  0  then  state  2  is  never  reached  in  the  Markov  chain 

with  transition  matrix  P^  violating  the  required  i rreducibi 1 i ty  property  of  P  . 

The  major  advantage  in  choosing  P  =ee'/m  is  that  the  irreducibi 1 i ty  of  each 

O  AWV 

P.  transition  matrix  is  guaranteed  at  each  step  of  the  procedure. 

The  next  stage  to  be  taken  with  the  procedures  proposed  in  this  paper  is  to 
carry  out  some  numerical  studies  and  to  compare  the  computational  speed  and 
accuracy  with  some  of  the  other  procedures  surveyed  in  the  initial  section. 
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